#!/usr/bin/python
# GM_review.py
'''Converts from cartesian galactocentric coordinates to earth-centric spherical coordinates.

Or that's the idea.  I know how to do it, in principal, just fine, but in practice couldn't figure a way that didn't involve four nested 'for' loops.  Ugly.  Doesn't work right now.'''

import numpy

def review(simpars, phypars, appearancexyz):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    # Relevant distances
#    x = numpy.array(xlocs, dtype=float)*PtScalexy
#    dx = numpy.subtract( x, R0 )
#    y = numpy.array(ylocs, dtype=float)*PtScalexy
#    z = numpy.array(zlocs, dtype=float)*PtScalez
#    dx2 = numpy.square(dx)
#    y2 = numpy.square(y)
#    z2 = numpy.square(z)
#    dR2 = numpy.add( dx2[:,numpy.newaxis], y2[numpy.newaxis,:] )
#    dR = numpy.sqrt(dR2)
#    d2 = numpy.add( dx2[:,numpy.newaxis], y2[numpy.newaxis,:] )
#    d2 = numpy.add( d2[:,:,numpy.newaxis], z2[numpy.newaxis,numpy.newaxis,:] )
#    d = numpy.sqrt(d2)
#    del x2, dx2, y2, z2, dR2, d2
    ds = dlocs*PtScaled
    las = llocs*PtScalela
    bas = blocs*PtScaleda
    zs = numpy.multiply(d[:,numpy.newaxis,numpy.newaxis], numpy.sin(b[numpy.newaxis,,numpy.newaxis:]))
    zindex = zs/PtScalez
    dz0 = zindex - numpy.array(zindex, dtype=int)
    dz1 = 1 - (zindex - numpy.array(zindex, dtype=int))
#    zindex = numpy.array(zindex, dtype=int)
    zindexlo = numpy.array(zindex, dtype=int)
    zindexhi = numpy.array(zindex, dtype=int) + 1
    dRs = numpy.multiply(d[:,numpy.newaxis], numpy.cos(b[numpy.newaxis,:]))
    xs = numpy.multiply( dRs[:,numpy.newaxis,:], numpy.cos( l[numpy.newaxis,:,numpy.newaxis] ) )
    xindex = xs/PtScalexy
    dx0 = xindex - numpy.array(xindex, dtype=int)
    dx1 = 1 - (xindex - numpy.array(xindex, dtype=int))
#    xindex = numpy.array(xindex, dtype=int)
    xindexlo = numpy.array(xindex, dtype=int)
    xindexhi = numpy.array(xindex, dtype=int) + 1
    ys = numpy.multiply( dRs[:,numpy.newaxis,:], numpy.sin( l[numpy.newaxis,:,numpy.newaxis] ) )
    yindex = ys/PtScalexy
    dy0 = yindex - numpy.array(yindex, dtype=int)
    dy1 = 1 - (yindex - numpy.array(yindex, dtype=int))
#    yindex = numpy.array(yindex, dtype=int)
    yindexlo = numpy.array(yindex, dtype=int)
    yindexhi = numpy.array(yindex, dtype=int) + 1
    del zs, ys, xs, dRs
    appearancedlb = []
    for i in appearancexyz:
        appearancedlb.append(numpy.ones(dres,lres,bres), dtype=float)
    emptiness = numpy.zeros(
    averageterms = numpy.array([0,0,0,0,0,0,0,0])
    weights = numpy.array([0,0,0,0,0,0,0,0])
    for i in range(0, len(appearancexyz)):
        averageterms[0] = appearancexyz[i][xindexlo, yindexlo, zindexlo]
        weights[0] = numpy.sqrt(dx0**2+dy0**2+dz0**2)
        averageterms[1] = appearancexyz[i][xindexlo, yindexlo, zindexhi]
        weights[1] = numpy.sqrt(dx0**2+dy0**2+dz1**2)
        averageterms[2] = appearancexyz[i][xindexlo, yindexhi, zindexlo]
        weights[2] = numpy.sqrt(dx0**2+dy1**2+dz0**2)
        averageterms[3] = appearancexyz[i][xindexhi, yindexlo, zindexlo]
        weights[3] = numpy.sqrt(dx1**2+dy0**2+dz0**2)
        averageterms[4] = appearancexyz[i][xindexhi, yindexhi, zindexlo]
        weights[4] = numpy.sqrt(dx1**2+dy1**2+dz0**2)
        averageterms[5] = appearancexyz[i][xindexhi, yindexlo, zindexhi]
        weights[5] = numpy.sqrt(dx1**2+dy0**2+dz1**2)
        averageterms[6] = appearancexyz[i][xindexlo, yindexhi, zindexhi]
        weights[6] = numpy.sqrt(dx0**2+dy1**2+dz1**2)
        averageterms[7] = appearancexyz[i][xindexhi, yindexhi, zindexhi]
        weights[7] = numpy.sqrt(dx1**2+dy1**2+dz1**2)
    return appearancedlb
